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CROSS REFERENCE TO RELATED APPLICATIONS 

[0001] This application is a Continuation-in-Part of U.S. Patent App. Ser. No. 

10/236,204 filed on September 6, 2002. App. Ser. No. 10/236,204 claims priority 
5 from United States Provisional App. Ser. No. 60/3 18,083 filed on September 7, 2001. 

App. Ser. No. 10/236,204 is a Continuation-in-Part of U.S. Patent App. Ser. No. 

09/580,863 filed on May 30, 2000, now U.S Patent 6,502,037. App, Ser. No. 

09/580,863 is a Continuation-in-Part of U.S. Patent App. Ser. No. 09/285,570 filed on 

April 2, 1999, now U.S. Patent 6,278,948; App. Ser. No. 09/399,218 filed on 
10 September 17, 1999, now U.S. Patent 6,424,918; and U.S. Patent App. Ser. No. 

09/405,850 filed on September 24, 1999, now U.S. Patent 6,430,507. App. Ser. No. 

09/399,218 filed on September 17, 1999, now U.S. Patent 6,424,918 and U.S. Patent 

App. Ser. No. 09/405,850 filed on September 24, 1999, now U.S. Patent 6,430,507 

are both Continuation-in-Part applications of U.S. Patent App. Ser. No. 09/285,570 
15 filed on April 2, 1999, now U.S. Patent 6,278,948. 

BACKGROUND OF THE INVENTION 
Field of the Invention 

[0002] The present invention pertains to processing potential fields data using vector 
20 and tensor data along with seismic data and more particularly to the modeling and 

inversion of gravity data which may also be combined with seismic data for 
hydrocarbon exploration and development. 

Related Prior Art 

25 [0003] Exploration for hydrocarbons in subsurface environments containing 

anomalous density variations has always presented problems for traditional seismic 
imaging techniques by concealing geologic structures beneath zones of anomalous 
density. Many methods for determining parameters related subsurface structures and 
the extent of the highly anomalous density zones exist. 
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[0004] United States patent number 4,987,561 , titled "Seismic Imaging of Steeply 
Dipping Geologic Interfaces, issued to David W. Bell, provides an excellent method 



for determining the side boundary of a highly anomalous density zone. This patent 
locates and identifies steeply dipping subsurfaces from seismic reflection data by first 
identifying select data which have characteristics indicating that seismic pulses have 
been reflected from both a substantially horizontal and a steeply dipping interface. 
These data are analyzed and processed to locate the steeply dipping interface. The 
processed data are displayed to illustrate the location and dip of the interface. This 
patent, while helping locate the boundaries, provides nothing to identify the 
subsurface formations on both sides of the boundary. 

[0005] There have also been methods for identifying subsurface formations beneath 
anomalous zones using only seismic data to create a model and processing the data to 
identify formations in light of the model. By further processing reflection seismic 
data, the original model is modified or adjusted to more closely approximate reality. 

[0006] An example of further processing seismic data to improve a model is United 
States patent number 4,964,103, titled "Three Dimensional Before Stack Depth 
Migration of Two Dimensional or Three Dimensional Data," issued to James H. 
Johnson. This patent provides a method of creating a three-dimensional model from 
two dimensional seismic data. This is done by providing a method of ray tracing to 
determine where to move trace segments prior to stacking. The trace segments are 
scaled to depth, binned, stacked and compared to the seismic model. The model can 
then be changed to match the depth trace segments that will be stacked better, moved 
closer to their correct three-dimensional position and will compare better to the 
model. This patent uses a computationally extensive seismic process to modify a 
seismic model that may be inaccurate. 

[0007] One source of geologic exploration data that has not been used extensively in 
the past is potential fields data, such as gravity and magnetic data, both vector and 
tensor data and using potential fields data in combination with seismic data to provide 
a more accurate depth model or to derive a velocity model. 



[0008] Gravity gradiometry has been in existence for many years although the more 
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sophisticated versions have been held as military secrets. The measurement of 
gravity became more usable in the late eighteen hundreds as measuring instruments 
with greater sesnsitivity were developed. Prior to this time, while gravity could be 
measured, the gravity gradient due to nearby objects could not be reliably measured. 

[00091 It has been known since the time of Sir Isaac Newton that bodies having mass 
exert a force on each other. The measurement of this force can identify large objects 
having a change in density even though the object is buried beneath the earth's surface 
or in other ways out of sight. 

[0010] Exploration for hydrocarbons in subsurface environments containing 
anomalous density variations such as salt formations, shale diapers and high pressure 
zones create havoc on seismic imaging techniques by concealing geologic structures 
beneath zones that are difficult to image. Often the imaging problem is due to 
velocity anomalies directly related to the anomalous densities of these structures. By 
utilizing gravity, magnetic and tensor gravity field measurements along with a robust 
inversion process, these anomalous density zones can be modeled. The spatial 
resolution obtained from this process is normally lower than that obtained from 
reflection seismic data. However, models obtained from gravity and magnetic data 
can provide a more accurate models for the seismic processing. Using the potential 
fields data models as an aid in seismic depth imaging processing greatly enhances the 
probability of mapping these concealed geologic structures beneath the zones of 
anomalous density. 

[0011] Zones of anomalous density may also be associated with zones of anomalous 
fluid pressure. Typically, while drilling an oil or gas well, the density of the drilling 
mud must be controlled so that its hydrostatic pressure is not less than the pore fluid 
pressure in any formation along the uncased borehole. Otherwise, formation fluid 
may flow into the wellbore, and cause a "kick." Kicks can lead to blowouts if the 
flow is not stopped before the formation fluid reaches the top of the well. If the fluid 
contains hydrocarbons, there is a serious risk of an explosion triggered by a spark. For 
this reason, wellbores are drilled with a slight excess of the borehole fluid pressure 
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over the formation fluid pressure. 



[0012] A large excess of the borehole fluid pressure over the formation fluid 
pressure, on the other hand, is also undesirable. Fractures in the borehole wall may 
result in loss of circulation of the drilling fluid, resulting in stuck drill strings, time 
delays and greater costs. Serious formation damage may also occur that can decrease 
the amount of recoverable minerals. 

[0013] Pressure prediction is done by estimating certain key parameters that include 
the overburden stress or confining stress, which is defined as the total lithostatic load 
on a rock volume at a given depth, and the effective stress, which is defined as the net 
load on the grain firamework of the rock at a given depth. These two relations are 
then used in the Terzaghi effective stress law to estimate the fluid or pore pressure. 
Terzaghi*s law states that: 

Pc = Pe -I- Pp where: 

(Pc) = the confining stress 

(Pe) = the stresses bom by the grains, and 

(Pp) = the stress bom by the fluid. 

[0014] Some workers treat a special case of Terzaghi's law where the confining stress 
is assumed to be the mean stress as opposed to the vertical confining stress. It should 
be acknowledged that this difference exists, but that it does not effect the 
embodiments of the present invention as they will pertain to estimating the total 
overburden load, which can then be converted to either vertical confining stress or 
mean stress based on the stress state assumptions that are made. A prior art method 
for estimating confining stress is to use a density log firom a nearby calibration well 
and integrate the density data to obtain the overburden load. This calibration is 
applied firom the mudline down to depths usually beyond the depth of sampling to 
predict the overburden away firom the calibration well. 

[0015] It has long been recognized that velocities of seismic waves through 
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sedimentary formations are a function of "effective stress," defined as the difference 
between the stress caused by the overburden and the pore fluid pressure. A number 
of methods have been used to measure the seismic velocities through underground 
formations and make an estimate of the formation fluid pressure from the measured 
velocities. Plumley (1980) and United States Patent 5,200,929 issued to Bowers, (the 
'929 patent) describe a method for estimating the pore fluid pressure at a specified 
location. The method also accounts for possible hysteresis effects due to unloading of 
the rock formation. The method utilized a pair of sonic velocity-effective stress 
relations. One relationship is for formations in which the current effective stress is 
the highest ever experienced. A second relationship is used when the effective stress 
has been reduced from the maximum effective stress experienced by the rock and 
hysteresis must be accounted for. 

10016] The '929 patent uses density data from nearby wells or from a geologically 
similar well to obtain the overburden stress. In most circumstances, the overburden 
stress may be adequately described by general compaction models in which the 
density increases with depth, giving rise to a corresponding relation for the relation 
between depth and overburden. In the absence of well control, determination of the 
overburden stress even within a sedimentary column is problematic. Furthermore, 
there are circumstances in which the model of a density that increases uniformly with 
depth is not valid. In such cases, the assumption of increasing density with depth is 
violated and a different approach to estimation of the overburden stress is needed. 

[0017] There are several types of situations that may arise wherein a model of 

density increasing with depth and compaction is not valid. In the first case, there is a 

region of abnormally high density in the subsurface, usually of magmatic origin. The 

region could consist of an extrusive or intrusive volcanic material having relative 

density of 2.8 or higher. When such a formation is present within a sedimentary 

section where the relative density is typically between 2.4 and 2.65, the result is an 

increase in the overburden stress underneath the formation over what would be 

determined by prior art calculations. On the other hand, a region of abnormally low 

density may occur from salt bodies (2.10) or shale diapirs. In such a case, the 
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overburden stress is abnormally low compared to what would be determined by prior 
art methods. In either case, even if the effective stress could be determined from 
seismic velocity measurements, a formation fluid pressure determination based on a 
prior art density model would be invalid. 

5 

[0018] Figure 1 is a seismic section illustrating an area of interest having highly 
anomalous density zones such as salt domes. A rounded interface can be seen as 
masking the formations below. For this data set, the lower boundary cannot be 
determined by normal seismic data processing methods. 

10 

[0019] There is a need for a method to image subterranean formations which are 
responsive to potential fields and non-potential fields data. There is a need for a 
method to combine these data types to extract more usefiil information than either 
data type has provided before. There is a need for a method that can more accurately 
15 determine the density of the subsurface in three dimensions away from and deeper 

than the limits of density from well control. The present invention satisfies this need. 

SUMMARY OF THE INVENTION 

[0020] The present invention is a method for determining a parameter of interest of a 
region of interest of the earth. At least one component of potential fields data is 

20 measured at a plurality of locations over a region of interest including a subterranean 

formation of interest. The potential field data are selected from magnetic data and 
gravity data. An initial geophysical model is determined for the region including the 
subterranean formation of interest. For the model, geophysical tensor data is updated 
using a forward model at a plurality of locations using a High Order Compact Finite 

25 Difference method and may include conjugate gradient methods. The method may be 

used to determine potential field measurements within a source region that is difficult 
to formulate using standard intregal equation approach where the Green's fiinction is 
evaluated. Green's functions evaluations are expensive and require significant 
memory to store, where as in this formulation only requires storing a sparse stencil. 

30 A difference between the estimated model value and the measured value of the 
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potential field measurements are determined at the plurality of locations. Based on 
this difference the geophysical model is updated. The model is iteratively updated 
and compared to the measured data until the differences reach an acceptable level and 
the parameter of interest has been determined. 



BRIEF DESCRIPTION OF THE DRAWINGS 



[0021] The patent or application file contains at least one drawing executed in color. 
Copies of this patent or patent application publication with color drawings(s) will be 
provided by the Office upon request and payment of the necessary fee. The present 
invention and its advantages will be better understood by referring to the following 
detailed description and the attached drawings in which: 

Figure 1 is a seismic section of an area having anomalous density zones such as a salt 
dome; 

Figure 2 A illustrates a 1 9 point finite difference stencil with matrix coefficients; 
Figure 2B illustrates the finite difference stencil of Figure 2A displayed as a cube; 
Figure 3 illustrates a comparison of an analytic solution with numerically calculated 
model data; 

Figure 4 illustrates a comparison of an integral formulation to a differential 
formulation; 

Figure 5A illustrates a realistic model with parameters which approximate a salt 
diapir; 

Figure SB illustrates the finite difference forward modeling results of a realistic 
model with parameters which approximate a salt diapir of Figure 5A; 
Figure 6 illustrates borehole gravity apparent density depth dependence; 
Figure 7 illustrates that the calculation of apparent density from borehole gravity 
using only the vertical component of the field results in an anomalously high density 
value; 

Figure 8 illustrates the contribution of each tensor component in the apparent density 
in borehole measurements; 

Figure 9 is a flow chart illustrating geophysical model development; 
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Figure 10 is flow chart of an embodiment of the invention using constrained 
optimization; 

Figure 11 is flow chart of an embodiment of the invention using Laplace's Equation 
and/or the Equivalent Source method; and 
5 Figure 12 is flow chart of an embodiment of the invention using a Poisson's Equation 

method solution. 

DESCRIPTION OF THE PREFERRED EMBODIMENT 

[0022] The present invention provides a 3D forward modeling method to compute 
gravity and tensor gravity data using Poisson's equation for airborne, surface and 
borehole applications in exploration problems. The forward modeling is based on a 
differential equation approach for solving the Poisson's (or Laplace's) equation using 
high order compact finite difference methods or finite element methods. 

[0023] The formulation of the present invention makes it possible to build inversion 
processes to predict or determine subsurface parameters, structures and areas of 
interest so as to minimize the difference between any or all of the model fields and 
their measured counterparts. Since this inversion process is so robust and applicable 
to so many varied field measurements, it is usefiil for exploration at both regional and 
prospect scales. 

[0024] United States Patents No. 6,278,948, No. 6,424,918, No. 6,430,507, and No. 
6,502,037 and U.S. Patent Application 10/236,204 having the same assignee and the 
25 contents of which are fiiUy incorporated here by reference, disclose methods and 

apparatus in which gravity and magnetics inversion techniques developed for use with 
gravity. Full Tensor Gradiometry (FTG) gravity and magnetics can be used as a driver 
for building earth models that may be used with seismic data processing. Some of the 
disclosures of the '948, '918, '507 and '037 patents and the '204 application are 
30 included here for convenience. 
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[0025] These previous inventions utilize very robust formulations for calculating 
vector and tensor gravity and magnetic fields due to a parameterized model such as 
salt formations embedded in a known or assumed distribution of sediments. The 
embedded model is composed of an upper boundary, usually mapped from seismic 
data, and a parameterized lower boundary to be determined by an inversion process. 
The parameterized lower boundary first uses gravity and/or magnetic data to predict 
parameters. Then, seismic data can be combined to provide a depth image using the 
predicted parameters as the driver. Further, a velocity model can be derived using the 
predicted parameters. In the alternative, the seismic data may be used as an 
additional constraint in running the inversion process again. This process of inversion 
followed by seismic imaging followed by another inversion and seismic imaging step 
may be repeated until the results of the potential fields inversion and the seismic 
imaging processes converge to a single answer or a point of diminishing returns is 
reached. Otherwise, if the process results begin to diverge indicating that there is not 
a unique solution, the process is discontinued. 

[0026] The inversion process of the present invention demonstrates its strength in 
several areas. In the area of background sediment and subsurface formation 
properties, the density can have any horizontal variability including discontinuities. 
The depth variation of these properties is assumed to be a polynomial up to order six 
for density. Thus, for practical purposes, there are no apparent restrictions on the 
variability of density especially considering the inherent resolution of gravity 
analyses. 

[0027] Computations based on parallel computing techniques are very fast for three- 
dimensional models. The inversion process, setup is ideal for using multiple 
processors. Thus, the elapsed time previously required for large three dimensional 
inversion processes can be greatly reduced. This is especially important for regional 
studies where large data sets are required. 

[0028] For data preparation, the forward modeling portion of this inversion method is 

also used to correct the gravity and gravity tensor data for complex bathymetric 
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effects. Bouguer corrections are made assuming a depth dependant density in order 
to obtain better input gravity and tensor data for the inversion process. This 
represents a substantial improvement over the standard use of a constant Bouguer 
density. 

[0029] Some of the major problems confronting the application of any modeling 
technique applied to potential field data are the removal of the regional field and, 
thus, the isolation of the field believed to be associated with the scale of the model 
being examined. By properly filtering the data between successive iterations in the 
inversion process, predicted regional fields can be obtained for gravity and magnetic 
data sets simultaneously, thus allowing convergence to a common model structure 
that explains both the band-limited gravity and magnetic fields associated with the 
anomalous body. The resulting regional fields then can be modeled to predict deep- 
seated structures beneath the area of immediate interest. 

[0030] Part of the inversion process that dramatically improves the convergence and 
efficiency of the algorithm is a predictive filtering procedure that reconstructs the 
regional field from the inversion itself. At each inversion step, the inversion 
estimates the anomalous body model that is required to fit the data and compares this 
model response to the observed field. The difference between the model and the 
observed field is treated as a "residual" or "error", the long wavelength component of 
this error field is calculated and attributed to the regional field that must be accounted 
for in the regional model. This long wavelength residual is used to reconstruct the 
regional model and this reconstructed regional model is compared to the long- 
wavelength regional component that is removed early in the preprocessing of the 
potential fields data to make sure that the signal used in the inversion was properly 
separated between signal related to the anomalous bodies and that related to the 
regional field. 

[0031] An important aspect of this invention is that the observation positions can be 

distributed arbitrarily, and even be positionied inside the source area, thus removing 

the need to have the data in a gridded form on the same datum. This is important 
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when simultaneously inverting different data sets and using borehole gravimeter data. 
Also, data sets from various contractors may have very different acquisition 
characteristics including geographic grids. Additionally, individual field observations 
may be weighted in the inversion process according to their uncertainty. 

[0032] Gravity data are generally acquired for the exploration of oil and gas or other 
natural resources like minerals and ore bodies. The different data types that are 
commonly acquired depend on the exploration objectives. The different types are: 
surface gravity data, tensor gravity data, airborne gravity data and borehole gravity 
data. The goal is to obtain subsurface density information or structural information 
from these different data sets by forward modeling and inversion. In principle, the 
density anomalies in the earth produce signals recorded in these data sets. Therefore 
one of the objectives is to compute the earth gravity response to subsurface anomalies 
using forward modeling with a known or estimated density model. The current 
invention of solving the Poisson's equation (using a high-order compact finite 
difference method) allows us to generate the gravity and tensor gravity responses in 
the entire medium as well as in the air with a single forward modeling run. This 
allows us to compute surface and borehole gravity responses concurrently without 
any additional complication in forward modeling. Prior art practice has been to 
compute the gravity and tensor gravity responses using an integral formulation that is 
much more complex, especially in the regions containing the source. This also 
requires more memory and computation time. The forward modeling approach 
method of the present invention provides the flexibility to invert surface, airborne and 
borehole gravity data simultaneously to recover density distribution in the subsurface. 

[0033] As is well known by practitioners in the art, the divergence theorem may be 
applied to the gravitational field 'g': 

V s 

With no attracting matter within the volume, V • g = 0 . This can be written as 
VVC/ = V^C/ = 0. 
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[0034] Prior art methods for modeling of tensor gravity data are based on solving the 
integral formulation for the different components of tensor gravity and a component 
for the vertical attraction. These methods can be used easily when the point of 
observation is outside the source to calculate gravity components on the surface. The 
methods involve the integration of each differential volume into which a model is 
discretized. 

[0035] If there is a particle of mass 'm' within the volume 'v' then the surface 
integral is different than zero and equal to: 



This gives Poisson's equation 

where p is the source density and / is a constant (the gravitational constant in the 
case of mass and gravitational potential). Poisson's equation may be written 
alternatively as, 

V^C/(jc, z) = -4;z7/7(x, 3;, z) . 

[0036] The forward modeling of Poisson's equation using the method of the present 
invention provides the flexibility to generate data in the entire medium including the 
source regions and thus allows inversion of surface and borehole gravity to recover 
density distribution in the subsurface. This is accomplished using the High-order 
compact finite difference (HOCFD) method. From a modeling standpoint this allows 
for generation of gravity and tensor gravity fields for the surface as well as for 
borehole data. 

[0037] Rewriting Poisson's equation as 



The term x^^ is the error from the approximation. Central differences are used to 
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calculate the operator in the discrete form using 

^ • 

Performing an expansion then: 



at 2!cbc' 3!ac' 4!6bc' 5!cbc' 



ay 

6!cbc* 



and 



(p{x -Sk)- 4>{x) —dK-\r 



dx 2\dx^ 3!cbc' 4\dx* 5\dx' 
Adding these equations and dividing by the cell thickness leaves 



6\dx^ 



2\dx^ ^ dx* 12 dx" 360 
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[0038] If all the terms in the Taylor series expansion are used, a high order 
approximation results. Common finite-difference methods eliminate the higher order 
terms leaving a second order approximation as 

^" 2!ax' ■ 
The error approximation term then is given by 



d*^ ^ a*^ ^ d'*^ 



•o(&r) 



ax* ay* dz\ 
where the value of the fourth order derivatives is given by 

a' fd^A d' ^ ^'-^ ^'-^^ 



dx' 



ay 



a;c^ 



/■ 



ay ay 
ay' az' j 



[0039] Substituting the error approximation into the original discretization results in 
20 the "High Order Compact Finite Difference" equation: 

o iZ 



This equation corresponds to a 19-point stencil as illustrated in Figure 2 A and Figure 

2B. In this equation, the right hand side is made up of the usual approximation 
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Sc^ +Sy^ -hSz^ and the higher order contribution — 



left hand side is made up of the source function fy^ and the second derivative of the 



[0040] The 19-point High Order Finite Difference Stencil "L" is illustrated in Figure 
2A and Figure 2B encompassing all the nodes of the mesh located on the three grid 
planes which intersect the center node of the cube, but not the comer points of the 
surrounding cube. An expanded view of the stencil with its matrix coefficients is 
illustrated in Figure 2A. The 5 solid vertices of the 9 point planes 201 and 205 of the 
stencil are combined with the 9 vertices on plane 203 to make up a 19 point finite 
difference stencil illustrated in Figure 2B. This is a 19 point operator for equi-spaced 
grid points. Stencils for other grids are known in the art and may be incorporated into 
the method of the present invention. For example, a 27 point operator formulation 
(for an o{h^) method) is described in *A High-Order Compact Formulation for the 
3D Foisson Equation* by W.F. Spotz and G.F. Carey, Numerical Methods for Partial 
Differential Equations, 12, 1996, pp. 235-243. 

[0041] Stencils for arbitrarily spaced grids are known in the art and may be 
incorporated into the method of the present invention. These arbitrarily spaced grids 
will allow for varying boundary conditions. 

[0042] Poisson's equation V^U ^ -4nyp may be expressed in matrix form as 



where L is the Finite Difference diagonal dominant matrix. The finite difference 
matrix can be stored in a compact sparse row format that requires only three vectors, 
therefore reducing memory significantly over other methods. The Finite Difference 
system of equations can be inverted using Conjugate Gradient methods and their 
variants that require one matrix vector multiplication. A conjugate gradient method 
implementation works quickly for compact sparse row matrices compared to other 



source function, —\Sc' 
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methods. 



[0043] The conjugate gradient method is a method for finding the nearest local 
minimum (the smallest value of a set, function, etc., within some local neighborhood) 
of a function of n variables which presupposes that the gradient of the function can be 
computed. A gradient in three dimension is, for example, the vector sum of the partial 
derivatives with respect to the three-coordinate variables y, and z of a scalar 
quantity whose value varies from point to point. The conjugate gradient method uses 
conjugate directions instead of the local gradient for going downhill. If the vicinity of 
the minimum has the shape of a long, narrow valley, the minimum is reached in far 
fewer steps than would be the case using other methods, for example the method of 
steepest descent. 

[0044] As an example of an embodiment of the present invention, a discretized 
operator L in the expression LU = f may be computed and examined for comparison 
using an analytic function for which the exact solution is known. Here 
U = sin;ccsin;zvsin;zz, 

/ = — sin;ccsin;zysin7E, for 
Q G [0,1]' . 

[0045] Figure 3 illustrates a comparison of the analytic solution for the observed 
data from the potential U on the left-hand side 301, 303 with the calculated data on 
the right 305, 307. The lower panels 303 and 307 illustrate the potential U at the 
center of the volume. 

[0046] Figure 4A illustrates an example of forward modeling the gravity fields using 

where the density change Ap = -OAgr/cm^ . The integral formulation results 
of the model are illustrated in Figure 4A with the vertical gravity , and gravity 
tensors , , . The integral formulation results may be favorably compared 
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with the results of the differential formulation as illustrated in Figure 4B with the 
vertical gravity , and gravity tensors g^^ ygyy^gzz- Figure 4C illustrates the gravity 
tensors g^ygy^yg^y the integral formulation and Figure 4D illustrates the gravity 
tensors g^^igyz^gxy for the differential formulation of the present invention. It is 
important also to note the boundary condition that the potential (LO goes to zero as 
distance goes to infinity. Figure 4E illustrates the convergence of the system of 
equation using conjugate gradient method. It took 10 min to solve on a SUN blade 
machine with Va million nodes. 

[0047] Figure 5A illustrates a realistic model with parameters which approximate a 
salt diapir with a density contrast of 0.4 grams per cubic centimeter, as shown by the 
change in density scale on the right. Figure SB illustrates the results of the method of 
the differential formulation of the present invention used to reconstruct the diapir 
with the vertical gravity g^ , and gravity tensors g,^,>gyy,g^ and g^ , , g ^ . 

[0048] As the method of the present invention can be used inside the volume 

containing the sources and inside the sources, this method has application to borehole 

gravity data. Prior art borehole gravity technology calculates an apparent density 

firom the observed vertical gravity as follows: 

g,(zo)-g,(zo+Az) 
• 

Figure 6 illustrates this Az dependence schematically. Apparent density is 

calculated as: 

pred ^Z}-Q 
r^a ^ zz ' 

rTty 

[0049] As illustrated in Figure 7, calculation of apparent density from borehole 
gravity using only the component 702 of the field results in an anomalously high 
density value AV. Combining the effect due to G^^ and as illustrated in panel 
703 results in a more accurate value for the density anomaly as shown in panel 701. 
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The contribution of each tensor component to the apparent density is illustrated in 
Figure 8. Panel 801 illustrates component G^^ normalized by G^^ + Gy^ -f G^^ , Panel 

803 illustrates component normalized by G^^ + -h G^ . Panel 805 illustrates 

component G^ normalized by G„ + G_^ + G„ . 

[0050] When the tensor fields are modeled with the differential form of the present 
invention we get, as a bonus, the capability of doing upward continuations as part of 
the calculations. The upward continued field is a result from computing the forward 
model. 

[0051] In a preferred embodiment of the present invention the forward modeling 
method for tensor gravity data is implemented by solving Poisson's equation using a 
4*^ order accuracy finite difference scheme. This enables accurate computation of the 
tensor gravity components. This implementation has the power of modeling surface 
tensor gravity data as well as borehole gravity data. Using horizontal tensor 
components of borehole gravity data contributes to a better prediction of the apparent 
density. Accurate modeling of borehole gravity data allows for more efficient 
monitoring of oilfields over time. The method adapts easily to compute upward and 
downward continued potential fields. 

[0052] Other variations are within the scope of the present disclosure. For example, 
the method is readily adaptable to be used within an inversion procedure that uses 
integral equations to calculate a forward model. 

[0053] Gravity, tensor gravity and borehole data are acquired in the exploration for 
oil, gas and minerals. These data provide a means to image the subsurface, based on 
density anomalies in the earth. The high order compact finite difference methodology 



using this work allows us to obtain accurate gravity (third order accuracy i.e. 0\h^ j 
and tensor gravity (second order accuracy, i.e. o{h^ )) in the entire medium. This 
methodology significantly improves the accuracy over the conventional nine point 
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finite difference approaches and therefore is highly advantageous for the modeling 
these field. The finite difference stencil used in a preferred embodiment is a 19 point 
stencil. The resulting system of equations is stored in a compressed row sparse format 
for use with the conjugate gradient inversion techniques as described earlier. 

[0054] The gravity (G^ , , ) and gravity tensor components 

(G„ , Gyy y , G^ , G^ , G^ ) are computed by obtaining the gravitational potential in 

the entire discretized medium using the above method. Subsequently the gravity and 
the tensor components are evaluated using forward and central difference formulas. 
The values of the components that are staggered with respect to the original mesh 
(coordinate layout) are interpolated to the initial nodal coordinates using a 3D 
interpolation method. The method of the present invention also provides for a 
formulation that can handle unequal cell sizes to compute the potential. Adapting the 
present High Order Compact formulation for use with unequal cell sizes has the 
advantage of imposing the Dirichlet boundary condition (i.e. where the potential goes 
to infinity) with fewer cells. 

[0055] Referring now to Figure 9, this flow chart illustrates a method for providing 
an initial model upon which derivation of subsurface formations (the geophysical 
model) can be based. At block 12 reflection seismic data is received. In general, 
reflection seismic data is more reliable than other forms of geologic exploration data 
in most situations because of the extensive work and research that has been devoted 
to its use and application. 

[0056] At block 14, the reflection seismic data is used to derive the top portion of a 
geologic model. As stated previously, reflection seismic data is very reliable, at least 
in this situation, for determining a model containing the top of an anomalous density 
zones. Geologic boundaries beneath the top of an anomalous zone or boundary are 
not as easily modeled. While reflection seismic data is generally more reliable than 
other forms of geological surveying, anomalous density zones, such as a salt dome, 
obstruct, divert or highly distort seismic reflection information below the boundary or 
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zone. In this type of area, reflection seismic surveying becomes unreliable, as shown 
in the example of Figure 1. In some cases, the reflection seismic data becomes so 
unreliable that no useful information below the salt dome can be obtained from 
reflection seismic exploration. 

[0057] At block 16 data pertaining to the determination of the lower boundary is 
received. This data may take the form any potential fields data, both vector and 
tensor. In the formulation used herein, the type of data is of no concem since any 
combination of the above mentioned data can be processed simultaneously. Although 
these types of data generally provide less resolution than reflection seismic data, in 
the case of an anomalous density zone, potential field data may provide the most 
reliable available data. 

[0058] The data received are used to formulate the limits of the lower boundary for 
the geologic model block 18. The actual lower boundary derivation is done by 
predicting parameters representing the lower boundary, and this predicting maybe 
done by a-priori knowledge or through an inversion process. 

[0059] Although various inversion techniques can be utilized to determine the 
parameters (coefficients) representing the Ipwer boundary, one preferred embodiment 
of the present invention involves the successive inversion of a single coefficient at a 
time until all coefficients are determined. The total number of coefficients is set a- 
priori at block 18 based upon the minimum wavelength (maximum firequency) desired 
in the lower boundary. Typical lower boundaries may contain as many as nine 
hundred coefficients for three-dimensional models, thirty in each horizontal direction. 
For example, if xi and X2 are the spatial limits of integration in the x direction, and j^/ 
and 72 are the spatial limits of integration in the y direction, and half cosine series are 
used, the number of terms required in the x and y directions respectively for a 
minimum wavelength ofX min are: 
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Thus the total number of coefficients representing the lower boundary is rix multiplied 



[0060] At block 20, the first coefficient representing a uniform lower boundary 
(longest wavelength component) is predicted. This coefficient is based on the 
received data and represents the best fit possible with limited information. 

[0061] At block 22 another coefficient is added through the inversion process. At 
block 24 all coefficients added thus far are readjusted through single parameter 
inversion in the order in which they were added to the lower boundary (firom the 
longest to the shortest wavelength coefficients). This process continues until all 
coefficients (rix times riy) have been predicted. At decision block 26 a determination 
is made as to whether all coefficients have been predicted. If they have not, the 
program returns to block 22 where another coefficient is added. 

[0062] If all coefficients have been predicted, the program proceeds to block 28 
where the lower boundary is displayed and subsequently sent on to the forward 
modeling of the present invention. The display may take any form currently in use in 
the art, such as a monitor display, a hard copy printout or be kept electronically on 
tape or disc. 

[0063] Initial geophysical model development and the updating of geophysical 
models based the differences between measured and estimated gravity and gravity 
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tensor data are extensively covered in the disclosures of the '948, '918, *507 and '037 
patents and the '204 application which have been fully incorporated herein by 
reference. Further updating geophysical models based on integration of potential 
fields and seismic data is covered in these disclosures as well. 

[0064] A number of embodiments of the '204 application are described therein 
which are applicable for combination with method of the present invention. In one 
embodiment, the predicted parameters are combined with seismic data to obtain a 
depth image and a velocity model for delineating formations or parameters of interest. 
In a second embodiment, imaging of the seismic data is performed using a lower 
boundary determined from the initial inversion. This seismic result is mapped to 
determine the position of the base of the anomalous body and is compared to the 
model parameters to obtain a difference between the two. The parameters are 
adjusted to provide a best fit. A new model is formed based on the new predicted 
parameters, and the potential fields modeling (using HOCFD in the present case) and 
seismic imaging modeling results are iterated until the solution converges. 

[0065] Where an anomalous body or subterranean formation of interest is a volcanic 
intrusive or extrusive body that has higher density than the surrounding sedimentary 
formations, the overburden stress will be higher than would be expected in a 
comparable depth of sedimentary rocks. The present invention is able to account for 
the effects of zones of anomalous density by determining the actual density, whereas 
prior art methods assume that density increases uniformly with depth or varies with 
depth based on limited well information. 

[0066] One embodiment of the present invention includes a modification of the very 

fast inversion process which produces models using vector and tensor potential 

methods (applicable for both gravity and magnetic data). This modification involves 

solving the inverse problem by using a constrained optimization approach, for 

example, a constrained nonlinear inversion method. The formulation allows 

imposing constraints on the solution as well as incorporating any geologic apriori 

information in the inversion algorithm. The computational speed is achieved using a 
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conjugate gradient solver and using approximate sensitivity when required. The 
algorithm has the flexibility to allow an interpreter to carry out hypothesis testing, 
based on his geologic intuition or any other information about the model derived from 
any sources, for example, geophysical or geological data. 

[0067] The method of this embodiment applies a constrained nonlinear inversion 
method using a Gauss-Newton approach. The model is parameterized using, for 
example, a boxcar basis function. Then the forward problem is solved using the 
method described above. In the inversion algorithm we invert for all model 
parameters simultaneously and optimally constrain parameters by imposing soft 
bounds on the solution. This bounded minimization problem is solved using an 
interior-point method. Interior-point methods are well known in linear programming. 
The interior-point algorithm guarantees that the solution will remain within bounds, 
i.e. the inverted base of salt will always be greater than (or equal to) the top of salt 
and less than any prescribed upper bound. This formulation has many model 
parameters compared to the data, which leads to an underdetermined problem. 
Therefore, to address the non-uniqueness issues, we minimize a model objective 
function subject to fitting the data. The choice of this model objective function is 
generally geologically driven. For example, we can impose smoothness criteria in 
some regions of the model and also allow sharp features to be built in other regions 
like faults that may be based on knowledge or geological objectives. In essence this 
model objective function can be used to impose a constraint on the model based on 
apriori knowledge. For example, a topological constraint may be placed on the base 
of the salt model (or any other parameter) based on apriori knowledge, and geological 
or geophysical information of the subsurface. 

[0068] Since we are solving for all of the model parameters in this approach, the 
system of equations may require a large matrix system be solved. Using a least- 
squares conjugate gradient solver and using approximate sensitivity when necessary 
achieves an increase in computation speed. The approximate sensitivity requires less 
storage and allows the matrix to be stored in a sparse form with very minimal loss in 

the information content due to approximation. The noise in the data is handled using 
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data weighting matrices. Appropriate stopping criteria for the convergence of the 
inversion process may be implemented and based on the assumption that the noise is 
random and has a Gaussian distribution. Other distributions and parameterization 
functions may be easily incorporated into this methodology. 

[0069] The constrained nonlinear inversion method using a Gauss-Newton approach 
leads naturally to an addressing of resolution issues for the inverted model. In the 
regions that are not constrained by the data, the model parameters from inversion may 
not be realistic. To avoid this problem, a reference model may be built into the 
inversion algorithm. Thus if the regions are insensitive to the data, the recovered 
model from inversion will reflect the reference model. This can be used to identify 
regions in the inverted model that are demanded by the data and thus provide a better 
interpretability of the results. 

[0070] The method of this embodiment is fiirther explained with reference to the 
flow chart of Figure 10. Gravity data 401, which may include gravity tensor data, are 
acquired for input to the method of the invention. An initial geophysical model 403 is 
determined or derived from apriori knowledge and/or other available data. A model 
potential field 405 is determined from the initial geophysical model 403. The initial 
geophysical model can be made using the forward modeling of the High Order 
Compact Finite Difference (HOCFD) method of the present invention. A difference 
is determined 407 between the acquired gravity data 401 and the gravity data 
potential fields' response due to the geophysical model 405. If the difference 407 is 
outside of a selected minimum range 409, an inversion process 410 is implemented, 
for example through applying a constrained nonlinear inversion using a Gauss- 
Newton approach as outlined in the '204 application. The model 405 is 
parameterized using appropriate basis functions; model parameters are inverted for 
simultaneously; soft bounds are imposed on the solution; the minimization problem 
may be solved with an interior-point method. The non-uniqueness of the 
underdetermined problem is addressed by minimizing one or more model objective 
functions, for example, by imposing geological constraints based on apriori 
knowledge or interpretations. In this manner, the estimated gravity potential fields' 
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response 405 is compared with measured gravity data 401 and iterated through 410 
until the model differences 407 are within selected limits 409. 

[0071] Two other embodiments of the present invention have application to address 
prior art data noise problems with the use of gravity tensor technology. One of the 
major problems with tensor methods is the typically high noise level associated with 
most tensor surveys. In fact, in many cases the noise may swamp the signals believed 
to be attributable to geologic structures. Before the tensor data can be used for 
exploration purposes, a process must be applied to try to extract the signal from the 
measured tensor data. 

[0072] Embodiments of the present invention for application to noise problems may 
incorporate one or both of two techniques that may be applied separately or in tandem 
to drastically reduce the level of noise in all five tensor channels. These two 
techniques are the 'Laplace's Equation' method and the 'equivalent source' method. 
Barely interpretable datasets can be converted to datasets that may be almost 
'geologic' in nature. 

[0073] The method based on a solution of Laplace's equation works very well for 
merging the low spatial frequencies measured with a standard gravimeter with the 
entire spatial frequency spectrum measured with ftiU gravity tensors measurements. 
The Laplace's equation method utilizes a solution on an arbitrary surface not 
intersecting the region containing the sources causing the gravity and tensor response. 
As an example, the method may use a 2D half-cosine transform-like ftinction for the 
X and y dependence, along with the exponential ftinction for the z dependence. Thus, 
the gravitation potential V(x,y,z), can be represented as a series of products of the two 
ftinctions with coefficients A^^ where m and n refer to the spatial wave numbers Kx, 
Ky and Kz in the x, y and z direction respectively. Kz becomes a ftinction of Kx and 
Ky when this series is made to be a solution to Laplace's equation for V(x,y,z). The 
utility of this form of V(x,y,z) is that it can be easily included in a linear-inverse 
scheme to investigate the true signal in a set of gravity and gravity tensor 
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measurements. Since the gravity field and tensor can be derived from the 
gravitational potential V(x,y,z), the 's can thus be estimated by fitting the 

observed data using various linear-inversion algorithms. This allows the components 
of all data channels that satisfy Laplace's equation simultaneously to be isolated from 
the noise. Noise may be defined as that part of the measurements that does not pass 
this test. The entire survey can be condensed down to a single set of coefficients . 

These coefficients can be used for producing arbitrary grids of all channels. In 
addition, various operations can be performed such as upward and downward 
continuation, higher-order derivative maps, etc. Furthermore, 3D modeling such as 
that outlined previously can be greatly simplified using output from this process. A 
case in point is that the modeling program can use one channel such as the second 
vertical derivative tensor (G^^), since it is directly calculated from V(x,y,z), and 
therefore contains information from all field measurements. This fact drastically 
reduces CPU time and memory requirements for doing 3D modeling. Hence, the 
method provides the opportunity of applying more sophisticated modeling techniques 
to larger survey areas while using fewer resources. One fiirther point is that this 
method can be applied to randomly distributed observation locations since no FFT 
operations are used and the whole process is model independent. 

[0074] The Laplace's Equation method embodiment is further explained with 
reference to the flow chart of Figure 11. Potential fields data 401 are acquired. An 
initial geophysical model 403 is determined or derived from apriori knowledge and/or 
other available data. A model potential field as a solution to Laplace's Equation 
method 405 is determined from the initial geophysical model 403, using for example, 
half-cosine function transforms in x and jv, and an z variable exponential function. A 
difference is determined 407 between the acquired potential fields data 401 and the 
potential fields response due to the geophysical model 405. If the difference 407 is 
outside of a selected minimum range 409, a forward modeling process 420 is 
implemented by applying a Laplace's Equation solution approach and/or an 
Equivalent Source method. In this manner, the estimated or modeled potential fields 
response 405 is compared with measured potential fields data 401 and iterated 410 
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until the model differences 407 are within selected limits 409. 

[0075] Another method that may be used for data noise reduction is based on 
equivalent-source modeling. Equivalent-source modeling deals with the formation of 
an infinitely thin layer of spatially varying mass between the source(s) believed to be 
responsible for the measured gravity and gravity tensor fields, and the observation 
positions. In the case of marine surveys this layer of surface mass may be chosen to 
be coincident with the bathymetric surface. The method of forward modeling has 
been outlined previously. The main difference is that the surface density is set as the 
product of a two-dimensional surface density distribution function, Ms(x',y'), and a 
delta function with the argument (z'-Z(x',y')), where Z(x',y') is the bathymetry. This 
results in an analytical integration of z', thus leaving a numerical 2D integration over 
x' and y' for the model field calculations. In a linear inversion program, Ms(x',y') 
can be determined that best models the measured tensors and gravity simultaneously. 
Ms(x',y') is parameterized analogously to the base of salt as previously disclosed. 
Just as for the Laplace equation solution method, the residual fields that cannot be 
modeled are considered noise, and this method also can utilize randomly distributed 
observation points. The Laplace equation solution method and the equivalent source 
modeling method may be used in tandem. 

[0076] The Equivalent Source method is fiirther explained with reference to the flow 
chart of Figure 11. Potential fields data 401 are acquired. An initial geophysical 
model 403 is determined or derived firom apriori knowledge and/or other available 
data. A model potential field as a solution to the Equivalent Source method 405 is 
determined from the initial geophysical model 403, using for example, surface 
density set as the product of a two-dimensional surface density distribution function, 
Ms(x%y'), and a delta function with the argument (z'-Z(x\y')), where Z(x',y') is the 
bathymetry. A difference is then determined 407 between the acquired potential 
fields data 401 and the potential fields response determined from the geophysical 
model 405. If the difference 407 is outside of a selected minimum range 409, a 
forward modeling process 420 is implemented by applying an Equivalent Source 
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method. In this manner, the estimated or modeled potential fields response 405 is 
compared with measured potential fields data 401 and iterated 410 until the model 
differences 407 are within selected limits 409. 

[0077] A preferred embodiment of the method of the present invention is further 
illustrated with reference to the flow chart of Figure 12. Potential fields data 401 are 
acquired. An initial geophysical model 403 is determined or derived from apriori 
knowledge and/or other available data using HOCFD. A model gravity data potential 
field as a solution to the Poisson equation 405 is determined from the initial 
geophysical model 403, using for example, data values derived a locations equivalent 
to the potential fields' data 401. A difference is then determined 407 between the 
acquired potential fields data 401 and the potential fields response determined from 
the geophysical model using HOCFD 405. If the difference 407 is outside of a 
selected minimum range 409, a forward modeling process 430 is implemented by 
appljang any inversion method wherein the modeled gravity and/or gravity tensor 
data are inverted to obtain values consistent with the structures of the initial 
geophysical model 401, or the structures and their density distributions may be altered 
as outlined in embodiments above. In this manner, the estimated or modeled 
potential fields response 405 is compared with measured potential fields data 401 and 
iterated 410 until the model differences 407 are within selected limits 409. 

[0078] While there has been illustrated and described particular embodiments of the 
present invention, it will be appreciated that numerous changes and modifications will 
occur to those skilled in the art, and it is intended in the appended claims to cover all 
those changes and modifications which fall within the true spirit and scope of the 
present invention. 
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